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Abstract — Fast surface integral equation (SIE) solvers seem to 
be ideal approaches for simulating 3-D nanophotonic devices, as 
these devices generate fields both in an interior channel and in the 
infinite exterior domain. However, many devices of interest, such 
as optical couplers, have channels that can not be terminated 
without generating reflections. Generating absorbers for these 
channels is a new problem for SIE methods, as the methods were 
initially developed for problems with finite surfaces. In this paper 
we show that the obvious approach for eliminating reflections, 
making the channel mildly conductive outside the domain of 
interest, is inaccurate. We describe a new method, in which the 
absorber has a gradually increasing surface conductivity; such 
an absorber can be easily incorporated in fast integral equation 
solvers. Numerical experiments from a surface-conductivity mod- 
ified FFT-accelerated PMCHW-based solver are correlated with 
analytic results, demonstrating that this new method is orders 
of magnitude more effective than a volume absorber, and that 
the smoothness of the surface conductivity function determines 
the performance of the absorber. In particular, we show that 
the magnitude of the transition reflection is proportional to 
1/L^''+^, where L is the absorber length and d is the order 
of the differentiability of the surface conductivity function. 

Index Terms — nanophotonics, surface conductive absorber, 
boundary element method, surface integral equation, reflections. 



I. Introduction 

In this paper we describe an absorber technique for terminat- 
ing optical waveguides with the surface integral equation (SIE) 
method, which otherwise has difficulties with waveguides and 
surfaces extending to infinity. In order to attenuate waves 
reflected from truncated waveguides, we append a region 
with surface absorption to the terminations, as diagrammed 
in Fig. [U The transition between the non-absorbing and ab- 
sorbing regions will generate reflections that can be minimized 
by making the transition as smooth as possible. We show how 
this smoothness can be achieved with the surface absorber 
by smoothly changing integral-equation boundary conditions. 
Numerical experiments demonstrate that the reflections of 
our method are orders of magnitude smaller than those of 
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Straightforward approaches, for instance, adding a volume 
absorptivity to waveguide interior In addition, we show the 
asymptotic power-law behavior of transition reflections as a 
function of the length of the surface absorber and demonstrate 
that the power law is determined by the smoothness of the 
transition [ 1 1. 

Unlike the finite-difference or the finite-element volume- 
discretization methods, SIE methods treat infinite homoge- 
neous regions (and some other cases) analytically via Green's 
functions, and therefore often require no artificial truncation 
of space. Because SIE methods only require surfaces to be 
discretized, they can be computationally efficient for prob- 
lems involving piecewise homogeneous media. A popular 
SIE method is the boundary-element method (BEM) [2J- 
|6|, particularly since the development of fast 0{NlogN) 
solvers iTll- lfTOl . However, a truncation difficulty arises with 
unbounded surfaces of infinitely extended channels common in 
photonics. Fig. [T] is a general photonic device schematic with 
input and output waveguide channels. In order to accurately 
simulate and characterize the device, such as calculating 
its scattering parameters, formally, the transmission channel 
must be extended to infinity, requiring infinite computational 
resources. A more realistic option is to truncate the domain 
with an absorber that does not generate reflections. 

The key challenge is to design an absorber that both has 
small reflections and is also easily incorporated into an SIE 
solver. The best-known absorber is a perfectly matched layer 
(PML) 1TT1-1|T3, but a PML coiTesponds to a continuously 
varying anisotropic absorbing medium, whereas SIE methods 
are designed for piecewise homogeneous media. A similar 
problem arises if one were to simply add some absorption 
within the waveguides — in order to minimize transition re- 
flection, the absorption would need to increase gradually from 
zero HI, corresponding again to inhomogeneous media. One 
could also use a volume integral equation (VIE) |16| or a 
hybrid finite-element method in the inhomogeneous absorbing 




Fig. I. Schematic diagram of a photonic device with input and output 
waveguide channels, which must be truncated in a surface integral equation 
method. 
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region, but then one would obtain numerical reflections from 
the discontinuous change in the discretization scheme from 
the SIE to the VIE. Moreover, it has been proposed that an 
integral-equation PML can be obtained by varying the Green's 
function instead of the media |17|, but a continuously varying 
Green's function greatly complicates panel integrations and 
fast solvers for SIE methods. 

With SIE methods, one could, in principle, implement a 
scalar waveguide absorption in a piecewise homogeneous 
fashion as a discontinuous increase in absorption, but this 
will obviously generate large reflections due to the discon- 
tinuity of the medium and also the numerical reflection due 
to the discretization of the interfaces. We demonstrate this 
with a finite rectangular waveguide, to which an absorber 
with constant volume electrical conductivity and magnetic 
conductivity is attached, as shown in Fig. |2] The longitu- 
dinal cross-section of this arrangement is shown in Fig. |3(a)| 
To achieve small reflections, the intrinsic impedance of the 
absorber is matched to that of the waveguide. Hence, the 
volume electrical conductivity (Te and the volume magnetic 
conductivity g^ should satisfy Gp/g„-i — ei/ ^i, where and 
Hi are the permittivity and permeability of both the waveguide 
and absorber media, respectively. We quantify the reflection 
by use of the standing wave ratio (SWR), the ratio of the 
maximum field magnitude to the minimum field magnitude 
in the standing- wave region, evaluated on the waveguide axis. 
From the SWR, a reflection coefficient is then readily obtained 
as in a conventional transmission line. We calculated the 
field reflection coefficients in this way for a range of Ge and 
Gm, and the smallest reflection coefficient obtained was 12%. 
This value, which is listed in the first column of Table |I] is 
unacceptable for many design applications. In particular, the 
design of tapers [IS] requires field reflection from terminations 
in the order of 10^"^ or smaller Evidently, a more sophisticated 
way of terminating waveguides is called for. 

In this paper, we examine an alternative approach to 
absorbers, adding electrical conductivity to the waveguide 
surface rather than to the volume, via a delta-function con- 
ductivity on the absorber surface, as shown in Fig. |3(b)| The 
absorber's interior medium remains the same as the waveg- 
uide's, thus eliminating the need to discretize the waveguide- 
absorber interface, shown as a solid line in Fig. |3(a)| and a 
dashed line in Fig. |3(b)| This surface-conductivity strategy 
permits an efficient surface-only discretization, but at the same 
time allows for a smoothly increasing surface conductivity, 
thereby reducing transition reflections. Specifically, surface 
conductivity is easily implemented in SIE methods as it 
corresponds to a jump discontinuity in the field boundary 
conditions at the absorber surface. Since the SIE explicitly 
discretizes the surface boundary, continuously varying the field 
boundary conditions is easily implemented. 

As a preview of results to be described below, in Fig. |4] 
we show the numerically computed complex magnitudes of 
the electric field along the x direction inside a rectangular 
waveguide with several different surface absorbers. The sur- 
face absorber is in the region where xq < x < {xq + L) 
in which Xq is the position of the interface and L is the 
absorber length. The surface electrical conductivity in this 
region is given by Gi,{x) — (7q( ^^^° )^, where d = 0,1,2 
for constant, linear and quadratic profiles. The constant gq 
is chosen so that the total attenuation over the length of the 



waveguide 



absorber 




Fig. 2. A discretized dielectric waveguide with an absorber attached. 



TABLE I 

The Standing wave ratio (SWR) and field reflection versus the 

CONDUCTIVITY DISTRIBUTION OF THE ABSORBER 



Absorber type 


Volume 


Surface 


Conductivity profile 


Constant 


Constant 


Linear 


Quadratic 


SWR 


1.2933 


1.1477 


1.0062 


1.0003 


Reflection r 


0.1279 


0.0688 


0.0031 


1.4998 X 10"* 



absorbing region matches that of the optimal volume absorber 
above. The approach for calculating the attenuation along the 
absorber will be explained in section IIV-AI In the complex 
magnitude plots of Fig. H] the peak-to-peak magnitudes of 
ripples are an indication of the amount of reflections. As is 
easily seen in Fig. IH there are substantial reflections when 
using a constant conductivity, smaller reflections when using 
a linearly increasing conductivity, and almost no reflections 
for a quadratically increasing conductivity. The magnitudes of 
the field reflection coefficients r are listed in Table |I] and show 
that the reflection coefficient for the quadratically varying 
surface conductivity is nearly one thousand times smaller than 
the reflection coefficient for the volume absorber The results 
for the constant, linearly and quadratically varying surface 
absorber verify the results in [1|, that the smoothness of any 
transition in the waveguide largely determines the resulting 
reflection. 

This paper is organized as follows. In section [III the BEM 
formulation incorporating surface conductivity is derived. In 
section |III1 the decay rate due to the surface conductivity 
is examined using both numerical experiments and calcula- 
tions using perturbation theory and Poynting's theorem. In 
section IIVI the asymptotic behavior of the transition reflection 
with respect to absorber length is presented. In section |V] 
we note the existence of radiation modes originating from 
the excitation source mismatch with the waveguide mode, and 
show that the radiation complicates the interpretation of certain 
numerical results, but does not affect the performance of 
the surface absorber. Details of the numerical implementation 
of the SIE solver are given in Appendix, specifically, the 
detailed matrix construction, acceleration and preconditioning 
techniques. 

ii. bem formulations with the surface 
Conductive Absorber 

In this section, we describe the 3-D BEM formulation for 
a waveguide truncated with a surface conductive absorber 
Fig. |3(b)| shows the x-z plane cross-section of an x-directed 
truncated rectangular waveguide, where the surface conductive 
absorber region is to the right of the dashed line. The permittiv- 
ity and permeability of the waveguide interior and the exterior 
media are denoted as e^, ji^ and e^, ^i, where the subscripts 
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(a) A waveguide with a volume absorber 




(b) A waveguide with a surface absorber 



Fig. 3. The 2-D longitudinal section of a waveguide with an absorber. 
The lengths of the waveguide and absorber are 20Xi and 10 A^. respectively, 
with Xi denoting the wavelength in the waveguide medium. The waveguide 
cross section size is 0.7211Ai X 0.7211Ai. The relative permittivities of the 
waveguide (silicon) and the external medium (air) are 11.9 and 1, respectively. 



e and i denote the exterior and the interior, respectively. The 
electrical surface conductivity ae (r) is subscripted with e as a 
reminder that only electrical conductivity is being considered, 
though the generalization of what follows to both electrical and 
magnetic conductivity could be considered. As is described 
in section Hill using only electrical conductivity can have a 
saturation phenomenon that can be avoided at the cost of using 
a longer absorber. The system is excited by a Gaussian beam 
propagating in +x direction. The Gaussian beam is generated 
by a dipole in a complex space lfT9l . where the real part of 
the dipole position is inside the waveguide, jXi from the left 
end. In this paper, the convention of the e^'^* time-harmonic 
mode is adopted. 

In SIE methods, for computing time-harmonic solutions, the 
unknowns are surface variables. In our case, we use surface 
electric and magnetic currents on both the interior, Jj and 
M,, and exterior, Je and Me, of every surface. The currents 
on surfaces with (Te = on the left side of Fig. |3(b)[ satisfy 
a simpler set of equations than the currents on surfaces where 
(Te 7^ 0, the right side of Fig. |3(b)| When appropriate, we 
distinguish between the Ce = and ^ currents with 
superscripts L and R, respectively. 

Invoking the equivalence principle i20l yields a relation 
between surface currents and fields 



— n X Ee = Me, 

h X He = Je, 

nx(E, +Einc) = M„ 

-h X (Hj +Hinc) = Ji, 



(1) 

(2) 
(3) 
(4) 



where n is an exterior-directed normal unit-vector, Ejnc, Hjnc, 
are the electric and magnetic fields of the Gaussian beam in a 
homogeneous space with material parameters equal to those of 
the waveguide interior, and Ee, He and E^, H^ are electric and 
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(c) quadratic surface conductivity along x direction 

Fig. 4. The complex magnitude of the electric field inside a waveguide and 
a surface absorber. The dashed line indicates the position of the waveguide- 
absorber interface. 
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magnetic fields due to the equivalent currents in the exterior 
and interior, respectively. 

On the surface of the waveguide, the continuity of the 
tangential components of the electric and magnetic fields 
yields the well-known PMCHW formulation 13, lE), 

nxE^(jf,Mf',jf,Mf) 

=n X [Ef (Jf , Mf , Jf , Mf ) + , (5) 

=n X [Hf(jf ,Mf, jf,Mf) + Hfj, (6) 

where E^(-) and H^(-) are integral operators described in the 
Appendix. From ([TJ-ll?]! and the tangential field continuity in 
the surface conductivity free region, the equivalent currents 
on the two sides of the waveguide surface are of the equal 
magnitude, but are opposite in direction. Specifically, 

Je = -Jf = J"^, (7) 

=-Mf = M^. (8) 

Thus, the unknown currents on the waveguide side are reduced 
to and M^. 

For the surfaces where de 7^ 0, a modified surface for- 
mulation is needed, one that incorporates the discontinuity 
due to the surface conductivity. When 7^ 0, the tangential 
electric field is still continuous across the absorber surface, 
and therefore 

nxEf (jf,Mf,jf,Mf) 

=h X [Ef (Jf , Mf , Jf , Mf ) + Eg,]. (9) 

The tangential magnetic field is not continuous as a sheet 
of surface electric current Ji„d = '^e^tLn is induced due 
to the electrical surface conductivity, thus creating a jump. 
Therefore, 

nx[Hf(jf,Mf,jf,Mf) 

-Hf (Jf , Mf , Jf , Mf ) - H«J = aeEf,„, (10) 

where E^,j — ~n x (n x E-'*^) is the tangential electric field 
on the absorber surface, and could choose the field on either 
side of the absorber surface according to the enforced equality 
in ©. 

As a result of ([T]i, (O and (|9]l, the interior and exterior 
magnetic currents can be represented with a single variable 

Mf = -Mf = M'^, (11) 

but the discontinuity of tangential magnetic field implies 
|Jfi 7^ |Jf I' ™d (01, (IHi, and ([TOll must be combined to 
generate a local equation 

Jf +Jf = aeEf,„. (12) 

Finally, using the integral operator relation between E^,-, and 
J and M, and substituting into (fTOl i. (fT2] i. yields 

n X [Hf(J^M^, jf,Mf) - Hf (Jf ,Mf , Jf ,Mf) 

-HL] =aeEf,,„(jf,Mf,jf,Mf), (13) 
Jf + Jf = ae[ Eft.n(jf , Mf , Jf , Mf ) + Ef (14) 

The unknowns, the surface currents J^, M^, J^, Jf 
and M^, can be determined by solving equations (|5]l, (|6]l 
on the left waveguide surfaces and (|9]l, ( fTsl l. (O on the 



Me 








J, M, 













Fig. 5. The 2-D longitudinal section of a waveguide with uniform sur- 
face conductivity. The waveguide length is lOAi and cross section size is 
0.7211Ai X 0.7211Ai. The relative permittivity of the waveguide and external 
medium are 11.9 and 1, respectively. 




Fig. 6. The complex magnitude of the electric field along x inside the 
waveguide in Fig. |5] with uniform surface conductivity. 



right absorber surfaces. The resulting linear system can be 
constructed and solved as described in Appendix. It is possible 
to simplify the formulation by eliminating one extra variable 
of electric currents on the absorber surface, though there are 
some accuracy issues ifSTl . 

III. The decay rate of the fields 

In this section the exponential decay rate of waves prop- 
agating through the surface absorber region is analyzed. We 
demonstrate the relation between decay rate and surface con- 
ductivity using the example of a single dielectric waveguide 
with uniformly distributed surface conductivity. The longitu- 
dinal cross-section is shown in Fig.|5] The behavior of interior 
fields generated by a Gaussian beam source were computed 
using a BEM method based on solving (|9]l, ( fTsT l and ( fT4] i. 

The plots in Fig. |6] show the complex magnitudes of the 
electric fields along the x axis inside the waveguide for two 
cases, (Je — 0.001 S and a^, = 0.002 S. As expected, the 
complex magnitude decays exponentially with distance from 
the source with a surface-conductivity dependent rate. Also, 
as can be seen, waves reflect back from the right end and 
presumably these reflections decay as they travel to the left. 

An approximation to the rate of exponential decay can be 
determined by fitting the field plots. The fitted decay rates 
for a range of Ce are shown in Fig. [7] and denoted with a 
dashed star curve. The decay rate does not monotonically 
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- * - exponential fit 

- ^ - first-order perturbation 

- - power conservation 
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Fig. 7. A comparison of three methods for computing the rate of field 
exponential decay along the propagation direction versus surface conductivity. 
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Fig. 8. Illustration of the approach using Poynting's theorem to calculate the 
decay rate of a waveguide with surface conductivity. The plot of the surface 
conductivity distribution cr^ix) along the longitudinal direction is aligned with 
the waveguide. 



increase with the surface conductivity. The curve shape can be 
explained as follows. When is small, the propagating wave 
is able to penetrate the lossy surface and is absorbed, with the 
absorption increasing with ag as expected. However, for large 
(Tg, the surface conductor itself becomes reflecting, forming 
essentially an enclosed metallic waveguide; as (Te oo the 
tangential electric field vanishes at the surface and therefore 
there is no absorption. The practical implications of this upper 
bound on effective values for de are limited, and are described 
in section HV] 

The following sections introduce two alternative approaches 
to calculate the decay rate from surface conductivities, to con- 
firm and further illustrate the numerical observations above. 

A. Decay rate calculation by perturbation theory 

In this section, a first-order closed-form decay rate formula, 
valid for low surface conductivity, is derived using perturbation 
theory. 

Assume the electric field E^°)(r) of the fundamental mode 
of a lossless rectangular waveguide is given, and the super- 
script (0) denotes the unperturbed quantity. When electrical 
surface conductivity cTg is put on the waveguide surface, it 
is equivalent to a perturbation of permittivity, denoted as 
Ae(r) = —'^6s{v), where 5s is the Dirac delta function 
across the waveguide surface. According to ll20l . ll22]| . the first- 
order variance of angular frequency due to the perturbation of 
permittivity Ae(r) is 



Aw 



(1) 



w /^A6(r)|EW(r)pdF 
2 /^e|E(0)(r)Pdy ' 



(15) 



where V is the whole volume domain and the superscript (1) 
denotes a first-order approximation. After applying the triple 
product rule to partial derivatives of the three interdependent 
variables lo, k, a^, we obtain a first-order change in propaga- 
tion constant k due to the frequency change in ( fTSl l. denoted as 
Afc(i) = where Vg is the group velocity, as Vg 

in which the subscript indicates ae is held fixed. Combining 
( fTSI l and equations of Ae(r) and Ak^^^ above, the integral in 



the numerator of (flST l is reduced to a surface integral of the 
tangential components of the electric field, therefore 



Afc(i) 



(16) 



2K,/^e|E(0)(r)|W 

where S is the surface of the waveguide. As expected, the 
perturbation in the propagation constant is imaginary, which 
in turn represents the decay rate a^^^ = —lm{Ak'^^^}. With a 
uniform cross section, the volume and surface integrals in ( fT6l ) 
can be further reduced to surface and line integrals on the cross 
section, respectively. The electric field before the perturbation 
E'^'^)(r), along with Vg, can be obtained numerically, for 
example, with a planewave method Ii23l . 

The decay rate calculated using the perturbation is plotted 
in Fig. |2] with a dashed diamond curve. Note that the curve 
overlaps with decay rates computed with other methods when 
the surface conductivity is small, and deviates for larger 
conductivity as should be expected given the first-order ap- 
proximation. 

B. Decay rate calculation using Poynting's theorem 

The perturbation method above predicts the decay rate when 
the surface conductivity is small. An alternative approach, 
based on Poynting's theorem, can be used to calculate the 
decay rate for the entire Ue range. 

Figure [8] shows a waveguide with surface conductivity and 
also plots the conductivity function af, {x) with x-axis aligned 
with the a;-axis of the waveguide. Since this approach requires 
integrating the fields of source-excited propagating modes in 
the exterior region, some inevitably excited modes, such as 
radiation modes that will be discussed in section [V] must be 
suppressed. For this reason, the surface conductivity starts with 
a large constant value (5 S). The large surface conductivity 
leads to the saturation as seen in Fig. |7] and therefore, the 
Gaussian-beam source excites metallic waveguide modes at the 
beginning, propagating in +x direction in the closed interior 
region. The surface conductivity is then smoothly reduced to 
a smaller value, with which the decay rate is to be calculated. 
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In this way, the metallic waveguide modes can smoothly 
change to the desired decaying dielectric waveguide modes 
with minimal radiation modes excited. 

The decaying propagating mode in the domain of interest, 
shown in Fig. |8] is assumed in the form 



E(r) 
H(r) 



Ho(y,z)e-"--^■''^ 



(17) 
(18) 



where /3 is the real propagation constant, and a is the unknown 
decay rate. The Poynting vector in frequency domain is given 
by S = -^ E X H*, and together with the assumed forms of 
E in ( 1171 ) and H in ( fTsT i, the derivative with respect to x is 
given by 

? = -2aS. (19) 

ax 

As illustrated in Fig. [8] apply Poynting's theorem in the closed 
volume V 

1 



RsJ S hdS^--J aelEtanpd^', (20) 

where S is the surface of the volume V, h is an exterior- 
directed normal unit-vector, and S" is the waveguide surface 
within V. In the limit as Ax — > 0, the closed integration 
surface S becomes a surface A, and one component of the 
integrand of the left side of (l20l i becomes ^ • x. Combining 
(fT9T l and this Aa; — )■ limit of ( |20] | yields a closed-form 
representation of the decay rate a 



Re /^(E X H*) • Ml 



2Re J^{E, X n*) ■ xdS 



(21) 



where L denotes the boundary of the surface A, and L' denotes 
the integral line on the waveguide surface within the surface 
A. 

The decay rate calculated using (1211 1 is plotted in Fig. |7] with 
a dashed circle curve. It shows good agreement with the decay 
rate computed using fitting for the entire de range, verifying 
that the surface conductivity is handled correctly by the BEM 
in accordance with Maxwell's equations. 

IV. Reflections and Asymptotic Convergence 

In section U we showed that a smoothly varying surface 
conductive absorber easily implemented in the SIE method is 
orders of magnitude more effective at eliminating reflections 
than a volume absorber of comparable length. And, since 
computational cost increases with the length of the absorber, 
it is worth examining the relationship between reflections and 
the absorber length. In order to do this, we first require a 
more careful classification of reflections and a more sensitive 
numerical measure of reflections than the standing wave ratio 
method. In the following subsections, we will define round- 
trip and transition reflections and present a measure of the 
power-law asymptotic convergence of transition reflections. 

A. Round-trip and transition reflections 

Reflections in a domain of interest can be divided into a 
round-trip reflection, Rr, and a transition reflection, Rt- The 
round-trip reflection is generated by waves entering into the 
absorber, propagating to the end without being completely 
absorbed, reflected off the end of the absorber, and eventually 



propagating back into the waveguide. The round-trip reflection 
coefficient is proportional to 



(22) 



where the coefficient C is determined by the reflection at the 
end of the absorber and L is the absorber length. A factor of 2 
in the exponent of (|22] | represents the effect of the round trip, 
and another factor of 2 indicates that the power is considered. 

The transition reflection is the reflection generated by the 
change in material properties at the waveguide-absorber in- 
terface. Therefore, smoother material change at the interface 
produces smaller transition reflection as predicted by coupled- 
mode theory (T). 

If the round-trip reflection Rr in (l22l i is held fixed, the 
decay rate a is inversely related to the absorber length L. In 
another words, for a longer absorber, a smaller decay rate a 
yields the same round-trip reflection. Since a is proportional 
to the surface conductivity CTc for small cTg, using a longer 
absorber implies a smaller transition reflection. It was further 
shown in HI that, given a fixed round-trip reflection, the 
transition reflection decreased as a power law with increasing 
absorber length L. The power-law exponent is determined 
by the order of the differentiability of the medium (conduc- 
tivity) function. Suppose a conductivity function along x is 
(Je{x) = aQu{ ^^j^° ), where the interface is at a; = Xq, and 
u{s) is a monomial function with order d in the < s < 1 
region 



u{s) 



< s < 1 
s < 



(23) 



With the surface conductivity function cre{x) and a fixed 
round-trip reflection, the asymptotic behavior of the transition 
reflection in terms of the length of the absorber is 



Rt{L) ^ O 



1 



1^2(1+2 



(24) 



The power-law behavior in ( |24] | indicates that, with a higher- 
order conductivity function, the transition reflection decreases 
faster with increasing the absorber length. It does not follow 
that d should be made arbitrarily large, however, there is a 
tradeoff in which increasing d eventually delays the onset L 
of the asymptotic regime in which (24) is valid ||T|. 



B. Asymptotic Convergence with L 

We present numerical results to verify the asymptotic power- 
law convergence of the transition reflection of the surface 
absorber. Because it is hard to explicitly measure the transition 
reflection in the integral equation method, instead, we measure 
alternative electric field expressions as in |1|. 

First, we define E(L) as the electric field at a fixed position 
in the waveguide when the length of the absorber is L, with 
unit \i. Thus E(L) includes the incident field, the round- 
trip reflection and the transition reflection. With a small fixed 
round-trip reflection, the difference of E(L + 1) and E(L) is 
the difference of the transition reflections, which in the limit 
of large L approaches to zero, so E(L + 1) — E(L) — > 0. 
Therefore, |E(i + 1) - E(i)|, and similarly |E(2L) - E(i)|, 
can be a measure of transition reflection, and specifically, 
derived from (l24l i. they are subject to the following asymptotic 



ZHANG et at: A NOVEL BEM USING SURFACE CONDUCTIVE ABSORBERS FOR FULL-WAVE ANALYSIS OF 3-D NANOPHOTONICS 



behavior 



|E(L + 1)-E(£)p 
|E(L)P 
|E(2L)-E(L)p 



O 
O 



1 



\ l,2d+2 



(25) 
(26) 



In the example of a rectangular waveguide attached with a 
surface absorber, the asymptotic convergence of dZSl ) is shown 
in Fig. |9(a)| in a log-log scale for constant, linear, quadratic, 
and cubic surface conductivity profiles {d — 0,1,2,3). The 
curve for d = slowly approaches to the L^^ curve. The 
other three curves align with the expected L^^,L^^,L^^^ 
curves respectively. The curves would eventually converge to 
a small quantity, which is the difference of the small round- 
trip reflections due to a phase difference. Fig. |9(b)| shows the 
asymptotic convergence of ( |26] | with the same waveguide ex- 
ample for constant, linear, quadratic, and cubic profile surface 
conductivities. Again, the curve for the constant conductivity 
profile converges slowly, and is expected to be aligned with 
the curve. The other three power-law convergence are 
the same as the predicted L~'^ , , . The agreement of 
the figures verifies the power-law behavior of the transition 
reflection of the surface absorber and further illustrates that a 
higher-order conductivity function leads to smaller transition 
reflections. 



V. Radiation in the surface absorber 

A surface absorber with a conductivity that rises smoothly 
with distance from the waveguide-absorber interface would 
be expected to have interior fields whose magnitude decays 
with distance at an accelerating exponential rate. Instead, the 
field magnitudes decay exponentially near the interface, but 
then switch to an O(^) decay, as shown in Fig. [10] The 
cause of this switch in decay rate is due to coupled radiation. 
The radiation is generated by the inevitable mismatch between 
the excitation source and waveguide modes. The situation is 
similar to a source radiating in a lossy half-space, in which the 
dominant field contribution is due to a lateral wave that decays 
algebraically Il24l . ||25l . Fig. [TOl shows the complex magnitude 
of the interior electric field in a waveguide attached with a long 
quadratic-profile surface absorber The waveguide system is 
excited by a dipole source and a Gaussian beam, respectively, 
located inside the waveguide. In the semilog plot, the dominant 
waveguide mode decays at an accelerating exponential rate, 
and then because the guided mode decays faster than the O(^) 
radiation, the radiation dominates at a certain distance from the 
waveguide-absorber interface. The coupled radiation results in 
an 0(i) floor, as is clearly shown in the inset with a log-log 
scale. Note that using a Gaussian beam source results in a 
lower floor than using a dipole source, because the Gaussian 
beam is more directional, and generates less radiation. 

The coupled 0(i) radiation does not affect the performance 
of the surface absorber, because the coupled radiation itself 
is several orders of magnitude smaller than the propagating 
modes in the waveguide, and little wiU be reflected. The 
asymptotic convergence of the transition reflections in sec- 
tion|lV]and the 10^* reflection coefficient in section U consis- 
tently showed the excellent performance of surface absorbers, 
obviously unhampered by the effects of radiation. 




(a) The asymptotic convergence of the transition reflection measured using 

|E(L-H)-E(L)p 



const 

- linear 

- quad 

- cubic 




10 10 10 

(b) The asymptotic convergence of the transition reflection measured using 

|E(2L)-E(L)|^ 

miw ■ 

Fig. 9. Asymptotic power-law convergence of the transition reflection with 
the length of the surface absorber. The length of the waveguide is lOXi, 
with \i denoting the wavelength in the waveguide medium. The waveguide 
cross section size is 0.7211Ai X 0.7211Ai. The relative permittivities of the 
waveguide (silicon) and the external medium (air) are 11.9 and 1, respectively. 



VI. Conclusion 

We presented a novel numerical technique, a surface con- 
ductive absorber that is easily combined with the boundary 
element method, to eliminate the reflections due to the trunca- 
tion of infinite channels. We illustrated this technique with a 
dielectric optical waveguide and described the modified BEM 
formulation needed to allow varying the surface conductivity. 
We further discussed the non-monotonically increasing decay 
rate with the surface conductivity and presented methods 
to calculate the decay rate using perturbation theory and 
Poynting's theorem. We demonstrated that the surface con- 
ductive absorber is orders of magnitude more effective than 
the volume conductive absorber, and showed an asymptotic 
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1/x 
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Gaussian beam 
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Fig. 10. The complex magnitude of the electric field inside a waveguide 
and a long surface absorber excited by a dipole source and a Gaussian 
beam, respectively. The lengths of the waveguide and the absorber are lOAi 
and 30 Ai, respectively. The surface conductivity on the absorber increases 
quadratically. 



power-law convergence of the transition reflection with respect 
to the length of the absorber to verify that the smoothness of 
conductivity function determines the transition reflection, also 
noted in 1 1 1 . 

The major advantages of the surface conductive absorber 
are: (1) the varying surface conductivity is easily implemented 
in BEM and can significantly reduce the transition reflections; 
and (2) the volume properties of the absorber are the same 
as the waveguide, so there is no interior cross section to dis- 
cretize, eliminating a potential source of numerical reflections. 

In this paper, the surface conductive absorber is demon- 
strated with a simple example, a rectangular waveguide. 
Demonstrating that the method is as effective on more general 
structures is the subject of future work. 
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VIII. Appendix 

A. Construction of a linear system 

From section HH the five equivalent currents, J^, on 
the waveguide surfaces, and Jf , Jf , on the absorber 
surfaces, are to be determined by solving the equations (|5]), 
(|6]l, ( fTSl l and ( fT4l ). The currents are approximated with the 
RWG basis function |4| on triangular- meshed surfaces. 



m 

M = ^M™X,„(r'), 



(27) 
(28) 



where Xm(r') is the RWG function on the mth triangle 
pair, and J„i, M„i are the corresponding coefficients for the 
electric and magnetic currents, respectively. These unknown 
coefficients of the five equivalent currents assemble a vector 
X to be solved for, specifically, 

"1^. (29) 



X = Af ^ J« Jf 



Electric and magnetic fields are represented using the 
mixed-potential integral equation (MPIE) Q for a low-order 
singularity, with integral operators L and if as in ||6l 



E*(J,M) = -Z,L*(J) + /f*(M), 
H*(J,M) = -X/(J) - li*(M), 



(30) 
(31) 



where Zi ~ \/Jm/q is the intrinsic impedance, the subscript 
I = e 01 i denotes the exterior or interior region, and the 
superscript t = L or R denotes the waveguide or absorber 
surfaces. The integral operators on the mth RWG function are 
given by 



Xm(r') + 72VV'-X™(r'; 



JC?(X„)=- / V X G,(r,r')X„(r')d5", 
Js' 



G,{r,r')dS' 

(32) 
(33) 



where S' is the surface of the mth triangle pair, ki — (jj^/Juei 
is the propagation constant in region I, and Gz(r,r') is the 
free-space Green's function in region I 



Gi{v,v') 



-jh\r- 



47r r - 



(34) 



where r and r' are target and source positions, respectively. 

We employ Galerkin method |26 | on the waveguide by using 
the RWG function as the testing function on target triangle 
pairs. The tested L, K operators on the nth target triangle 
pair due to the mth source triangle pair become 

4„„.(X„0 = /x„(r).i*(X,„)d5, (35) 
Js 

/CL™(X™) = f X^{r)-Kf{X^)dS, (36) 
Js 

where S is the surface of the nth target triangle pair Substi- 
tuting the tested field operators into equations (|5]l, © yields a 
matrix All due to the currents J^, M^, and a matrix Alt? 
due to the currents Jf , Jf , M-'^. 

On the absorber surfaces, the term cre(r)E^j^j^ in ST3[ and 
(fl4l i requires another testing procedure in order to incorporate 
the surface conductivity 

-Cf^^CX™) - /ae(r)X„(r).Lf(X„)d5, (37) 
Js 

^fr™(X™) = / ae(r)X„(r).Af (X„0d5. (38) 
Js 

Similarly, substituting the above four tested integral operators 
into (|9]l, ( flJl l and (fT4l i generate a matrix Arl due to the 
currents J^, M^, and a matrix Afiji due to the currents , 
Jf , M«. 

Assembling the four matrices according to (|5]l, (|9]l, (fTST l 
and (fT4l l yields a dense linear system 



Ax = b, 



(39) 
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where 



All Alr 

AflL -^RR 



(40) 



The right-hand-side vector b contains tested incident electric 
and magnetic fields 



in which the nth entry of each subvector is given by 

bkn = / X„(r).Ef„,(r)d5, 
Js 

f X„(r).HLe(r)d5, 
Js 



"Ea,n 



[ ae(r)X„(r).EL(r)d5. 
Js 



B. Acceleration and Preconditioning with FFT 



(41) 

(42) 
(43) 
(44) 



The linear system (|39] l can be solved with iterative al- 
gorithms, for instance, GMRES for this non-symmetrically 
dense system. In each iteration, the matrix-vector product takes 
0{N'^) time, where N is the number of unknowns. Moreover, 
to explicitly store the matrix A is expensive, requiring 0{N'^) 
memory. In fact, there have been many well-developed fast 
algorithms to reduce the costs of the integral equation solvers 
Q-lTOl. In this paper, we use a straightforward and easily- 
implemented FFT-based fast algorithm to accelerate the SIE 
method on periodic guided structures. 

As shown in Fig. (TT] the waveguide is discretized into a 
periodically repeating set of the RWG triangle pairs. Due to the 
mesh periodicity and the space invariance of the operators ( l35T l, 
(|36] |. the matrices All and Alr are block Toeplitz, requiring 
explicit calculation and storage of only a block row and a 
block column, reducing memory to approximately 0{N). A 
Toeplitz matrix can be embedded in a circulant matrix, and the 
circulant matrix-vector product can be computed with the FFT 
|l27|, |28|. In this way, the computational costs are reduced 
approximately to O(A^logA^). 

Because the surface conductivity CTc must vary with distance 
from the waveguide-absorber interface, the tested potential 
operators ( |37] |. (l38T l are not space invariant. Therefore, ac- 
celerating the (Tg parameterized matrices Arl and Arr 
is not as straightforward as All and Alr. Typically, the 
integral of dSTl l. (|38] | is calculated numerically using Gauss 
quadrature [29], summing up the tested potentials at quadrature 
(target) points with Gauss weights. The space invariance of 
the potential operators (|32] |. (|33] | and periodicity of the mesh 
allows assembling a matrix of the potentials at target points by 
explicitly calculating only a block row and block column. Then 
the potentials at target points are summed after testing and 
multiplications with Gauss weights and surface conductivity, 
and eventually stamped into matrices A^l and Arr. 

Another great advantage of working with a Toeplitz or a 
block Toeplitz matrix is the existence of a highly efficient 
preconditioner |30|-|32|. A circulant matrix is approximated 
from the Toeplitz matrix, and then can be easily inverted with 
the FFT. We use this method to calculate a preconditioner 
for All, and use the block-diagonal preconditioner IfTOl for 




Fig. 11. A discretized waveguide with a periodic unit. 
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